clear; close all; clc
 addpath('../LibForSerial/LibrarySLP','../LibForSerial/LibraryLeSage/','../LibForSerial/LibVAR');
%----------- Load shocks using two identification schemes ----------------%
 NumberConsecutiveshocks = 2; % 1st previous shocks that is positive, and the 2nd one is positive too 

  shockchoice = 1; % 1 means BdBG ReStud shock (from VAR system with IP)
% shockchoice = 2; % 2 means BdBG ReStud shock (from Predictive regressions)
if     shockchoice==1
cd ./DataShocks/
    LoadShocksMandQpositiveRVsp500FromVAR;sample=[ 1983 6; 2019  12];
elseif shockchoice==2
cd ./DataShocks/
    LoadShocksMandQpositiveRVsp500;       sample=[ 1983 3; 2019  12];
%   uRV_from_PredReg = uRV; clearvars -except uRV_from_PredReg 
%   save uRV_from_PredReg uRV_from_PredReg
%   STOP
end
  load uRV_from_PredReg
fprintf('Correlation between uRV from predictive vs VAR is %3.3f \n', corr(uRV_from_PredReg(end-size(uRV,1)+1:end),uRV)); 
cd ../

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data = xlsread('../Dataset/VARdataBBupdatedMonthly.xlsx','MainData');
startsample=find(data(:,1)==sample(1,1)*100+sample(1,2));
endsample  =find(data(:,1)==sample(2,1)*100+sample(2,2));

data = data(startsample:endsample,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data(:,[3:10]) = log(data(:,[3:10]))*100;
%--------------------------------------%
clear startsample endsample 

T = size(data,1);

%P  =  4; % From BdBG "All VARs below use four lags, as suggested by the Akaike information criterion for the main specification."
 P  =  6; % including in the LP regression P lags of all variables in the system (for consistency througout)

 
  plottype = 'Output';      ii=3; %
% plottype = 'riskfr';      ii=1; %              - Figure B.1, Panel (c-d)
                           %ii=2; % for Wu-Xia shadow rate
% plottype = 'Inflat';      ii=9; % PCE          - Figure B.2, Panel (c-d)
% plottype = 'SP500' ;      ii=7; % real SP500   - Figure B.3, Panel (c-d)
%% IRF 
h1 = 0 ;  %start LP at h=0 or h=1 (h=1 if impose no contemporaneous impact)
H  = 12*3; %#horizon of IR

% We are interested in the estimation of the dynamic multiplier of y_t+h with respect to a change in x_t for h ranging from 0 to H
%------------ IRF of GDP growth to a shock variable ----------------------%
      y  = data(:,ii);    %endogenous variable
%-------------------------------------------------------------------------%
  Remove_Large_Shocks = 0; % Set it to 1 for Online Figure B.4 Panel (c)-(d)
    x  = uRV;      
    IndSeqShocks = IndPreviousNShocksRV;
    IndLargeShocks = uRV>2;

LargeShocksToBeRemoved = IndSeqShocks & IndLargeShocks;  
% [x IndSeqShocks IndLargeShocks]
if Remove_Large_Shocks; IndSeqShocks = IndSeqShocks - LargeShocksToBeRemoved; end
fprintf('Number of events with %d consecutive shocks is %d \n',NumberConsecutiveshocks,sum(IndSeqShocks));

%-------------------------------------------------------------------------%
xs = x.*IndSeqShocks; %shock variable times states
%-------------------------------------------------------------------------%
  w  = [lagmatrix( data(:,unique([ii 1 10  ])), 1:P )                    ]; %control for employment as in BdBG ReStud VAR; otherwise IP response to 1 shock is positive!
if ii == 9
  w  = [lagmatrix( data(:,unique([ii 1 10 6])), 1:P )];
end

w( ~isfinite(w) ) = 0;
OKforReg = isfinite(x) & isfinite(y) & isfinite(xs);
for cc=1:size(w,2)
    OKforReg = OKforReg & isfinite(w(:,cc)) ;
end  
  x = x(OKforReg,:); 
  y = y(OKforReg,:); 
  w = w(OKforReg,:);
  xs=xs(OKforReg,:);
  
% This shrinks the estimated dynamic multiplier to a polynomial of order r - 1
r = 3; %(r-1)=order of the limit polynomial (so r=3 implies the IR is shrunk towards a quadratic polynomial)
lambda =  20; %value of penalty
 if ii==8 || ii==9;  lambda = 10000; end                                      
 if ii==6 || ii==7;  lambda = 10000; end  
%-------------------------------------------------------------------------%
slp_sd  = locprojStateDep(y,x                  ,xs,w,h1,H,'smooth',r,lambda); %IR from Smooth Local Projection
 lp_sd  = locprojStateDep(y,x                  ,xs,w,h1,H,'reg');             %IR from        Local Projection

%% Figure 2, Panels (c) and (d)  State-dependent IR to a RVAR shock
figure(2)
%-------------------------------------------------------------------------%
subplot(1,2,1)
  plot(       0:H, (slp_sd.IR+slp_sd.IRstate), 'b--', 0:H, (slp_sd.IR),'k--', 'LineWidth', 1.5 ); 
  title(['IR_{slp} when there are ' num2str(NumberConsecutiveshocks) ' Consecutive Positive Shocks']);  legend(         's_t=+1 (Previous shock(s)=1)','s_t=0','Location','Best')
hold on; plot(0:H , zeros(H+1,1) , '-k' , 'LineWidth' , 2 ); grid; xlim([0 H]);                        %legend('Linear','s_t=+1 (Previous shock(s)=1)','s_t=0','Location','Best')
if ii==2 
    ylim([-2.8 1.6]);set(gca,'YTick',[-2.8:0.5:0.6],'FontSize',12);ylabel('Percent','FontSize',12);
end
xlim([0      H]);set(gca,'XTick',[0 3:3:H],'FontSize',12);

subplot(1,2,2)
grpyat=[(0:1:H)', slp_sd.Interval_up_state; (H:-1:0)' slp_sd.Interval_down_state(end:-1:1)]; patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]); hold on
plot(0  :H , zeros(H+1,1), 'k-'); hold on 
plot(    0:1:H,   slp_sd.IRstate, 'b--','LineWidth', 1.5); hold on
title('State multiplier of State-dependent IR to uncertainty shock')
grid;
xlim([0      H]);set(gca,'XTick',[0 3:3:H ],'FontSize',12);